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Abstract 

In this paper, we reexamine the vahdity of using time quantified Monte Carlo (TQMC) method 
[Phys. Rev. Lett. 84, 163 (2000); Phys. Rev. Lett. 96, 067208 (2006)] in simulating the stochastic 
dynamics of interacting magnetic nanoparticles. The Fokker-Planck coefficients corresponding to 
both TQMC and Langevin dynamical equation (Landau-Lifshitz-Gilbert, LLC) are derived and 
compared in the presence of interparticle interactions. The time quantification factor is obtained 
and justified. Numerical verification is shown by using TQMC and Langevin methods in analyzing 
spin-wave dispersion in a linear array of magnetic nanoparticles. 

PACS numbers: 75.40. Gb, 75.40.Mg, 75.50.Tt 
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I. INTRODUCTION 



The TQMC method is found to be a powerful simulation technique in modeling mag- 




It is found that 



netization reversal dynamics of magnetic nanoparticles 
simulation with the TQMC method is considerably more efficient than the conventional 
method of modeling magnetization dynamics based on time-step integration of the stochas- 
tic LLG equation, especially in the case of high damping limits |3| . The attraction of TQMC 
also lies on the fact that it establishes an analytical connection between the two stochas- 
tic simulation schemes, Monte Carlo (MC) and Langevin dynamics, which were previously 
thought to have different theoretical bases. Such analytical connection provides alternative 
techniques to both stochastic models, e.g. solving a stochastic differential equation using 
advanced Monte Carlo techniques to calculate the long-time reversal jj, ^. 

The validity of using TpMC to simulate an isolated single domain particle is first demon- 
strated by Nowak et al. \l\ and later rigorously proved by us in Ref. jsl using the Fokker- 
Planck equation as a bridge between MC and Langevin methods. In the case of interacting 
spin arrays, the validity of TQMC has not been analytically proved although it has been 



numerically shown 0, |5[ . It thus comes the necessity to establish the proof for the case of 
interacting spin systems, since in practical applications the discrete spins or moments (in the 
form of magnetic nanoparticles) are usually closely packed together, and hence are strongly 
coupled to one another. It is also important to show explicitly whether the analytical equiv- 
alence between the TQMC and the stochastic LLG equation and the time quantification 
factor, are dependent in any way on the nature (e.g. magnetostatic or exchange) or strength 
of the coupling interactions. 

In this paper, we jjrovide a rigorous proof for this case, based on the technique presented 
in our earlier works |5]. We further demonstrate the generality of the TQMC method, by 
implementing it in two different contexts, i.e. in time-evolution and reversal studies of a 
square array of spins, and in analyzing the spin wave dispersion in a linear spin chain. 

II. MODEL 

The physical model under consideration is a spin array (which represent an array of 
magnetic nanoparticles), whose spin configuration is represented as {s} = {■ ■ ■ ,Sj,---}, 
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where s = M/M^ is a normalized unit vector representing the magnetic moment of the spin 
and i refers to the i^^ spin in the vector hst of length N. 

The micromagnetic dynamics of the spin array is traditionally described by the Landau- 
Lifshitz-Gilbert (LLG) equation: 

is} = --^{s} X ({h} + a.{s}x {h}) (1) 
at 1 + 

where a and 70 are the damping constant and the gyromagnetic constant respectively, {h} = 
2|^-yV{s}-E is the effective field which is normalized with respect to the anisotropy field 
Hk = 2Ku/ fio^s, where is the anisotropy constant and /io is the magnetic permeability. 
E = E{{s}) is the total energy of the system which consists of the typical contributions 
in a micromagnetic system, e.g., Zeeman term, anisotropy term, magnetostatic term and 
exchange coupling term. The operator {s} x {h} = {•■-, Sj x hj, ■■■} is understood. To 
represent the thermal fluctuation, white noise-like stochastic thermal fields are added to the 
effective field according to Brown [3]. 

Alternatively, Random walk Monte Carlo (MC) algorithm can also be used in simulating 
the magnetization reversal dynamics ^1]. At each Monte Carlo step, one of the N spin 
sites is randomly selected to undergo a trial move, in which a random displacement lying 
within a sphere of radius R {R <^ 1) is added into the original magnetic moment and 
the resulting vector is then renormalized. The magnetic moment changes according to a 
heat bath acceptance rate as A{AE) = 1/(1 + exp{i3AE)). Here AE is the energy change 
within the random walk step and (3 = {ksT)'^, ks is the Boltzmann constant and T is the 
temperature in Kelvin. 

III. FOKKER-PLANCK EQUATIONS 

To link the MC scheme with the stochastic LLG equation, we shall derive the respective 
Fokker-Planck (FP) coefficients corresponding to the LLG equation and the random walk 
MC j^l. The general Fokker-Planck equation (FPE) for a spin array in a spherical coordinates 
is given as 

''.pm, {V}. = - E I- (^'. ■n-Y.ir (A,. ■p) + \Y.i^ {b,,. ■ p) 



19^ 1 (9^ 



where drift terms and diffusion terms B^y (x = {9i,(pi}, y = {9j,fj}) are defined as the 



ensemble mean of an infinitesimal change of x and y with respect to time 



By giving the 



detailed derivation in the appendix, we obtained the Fokker-Planck coefficients for LLG: 
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as well as for the TQMC: 



20 V 99, 
B^f, = N-^^ . S., (4) 

1 R2 
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where in Eq. /i' = ^„v-m7(iW) ' ^' = ^V«, = /^V/?- 



IV. MAPPING MC TO LLG 



In the high damping limit where the damping constant a is large, so that g' = h'/a 0, 
a term-wise equivalence can be established between the FPE coefficients in Eqs. and 
(jH), corresponding to the LLG and MC methods, if: 

Eq. (jni), in which Atmc is calibrated in MCS/site (one Monte Carlo step for each site 
on the average), is the time quantification factor for the TQMC method in interacting spin 
arrays. The time quantification factor is found to be the same as the one in Ref. [l| for 
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an isolated single particle case, and is thus consistent with previous numerical convergence 
observed in Refs. 0,0]. 

For the low damping limit where precessional motion becomes significant, one may wish 
to use the precessional (hybrid) Metropolis Monte Carlo algorithm ^5]. We confirm that, 
by using the same derivation techniques, one is able to prove the validity of including the 
precessional move in the MC algorithm in simulating the micromagnetic properties of an 
interacting spin array. 



V. RESULTS AND DISCUSSION 

The equivalence between the MC method and LLC, which is expressed by Eq. (0), pro- 
vides the theoretical justification for the use of MC method as an alternative to the LLG 
equation in micromagnetic studies. The equivalence which has been established is very gen- 
eral because no explicit form of the Hamiltonian is used in the derivation. This implies that 
the validity of the equivalence is independent of many physical and simulation parameters. 
For illustration, we test the validity of the TQMC method for a simple 10 x 10 spin array 
which is subject to a varying exchange coupling strength J . As shown in Fig. (^, the time 
evolution behavior of the (asymmetric) magnetization reversal is simulated for different val- 
ues of J . We find good convergence between the simulated results from both LLG and MC 
schemes, even when the switching mechanism of the spin array changes from the indepen- 
dent reversal (small J) to the nucleation-driven reversal (large J) . We also confirm that the 
mapping between MC and LLG time steps as expressed in Eq. (0), is also independent of 
other simulation and physical parameters, e.g. the chosen boundary condition (periodic / 
free), the lattice size, and the nature of the coupling (magnetostatic / exchange). 

Next, we show that the equivalence between the MC method and LLG enables the MC 
method to be utilized in most of the situations where LLG applies, and beyond the above 
time-evolution simulation. As an example, we consider the dispersion relation for the pri- 
mary spin wave mode of a one-dimensional spin chain. This example is chosen because it 
tests the capability of precessional TQMC method to simulate both spatial and time cor- 
relation of the spin-wave dynamics. By comparison, conventional MC methods are more 
suited for equilibrium or steady-state studies rather than time correlation dynamics. 
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FIG. 1: (color online) The time evolution behavior of the magnetization reversal in a spin array 
system. The following simulation parameters are assumed: lattice size of 10 x 10, periodical 
boundary condition, thermal condition KuV/ksT = 25, damping constant a = 1.0 and external 
field h = 0.5 applied on an angle 9 = it/4: with respect to the easy axes. The exchange coupling 
strength J is the adjustable variable. To guarantee the simulation accuracy, the time interval At 
for the LLG integration changes with J as At = 0.01/(1 + /i + J/KuV) |9|, while the trial move 
step size R in the MC simulation is chosen to reflect the At in one MCS. Error bars are smaller 
than the symbol size. 

The Hamiltonian of the spin chain system is set to be: 

H = ^ I - J ^ s, ■ % - K^V{s, ■ Kf - fioM.V ■ Si ■ Hext J (6) 

i y j<i 1 

where {i} represents the neighboring spins of the spin, J is the coupling strength, Hext 
is the apphed field and k„ refers to the unit vector along the easy axis. Magnetostatic 
coupling was not included in this test. The dispersion relation for the one- dimensional spin 
wave mode has been theoretically studied [3| and is given by: 



^(A:) = ^^^[1 + ^cxt + 4(J/2i^„\/) sin2(fca/2)] (7) 
1 + 



where /lext = -f^ext/-f^fc and a is the lattice constant. The calculations were done using the 
computational techniques of Refs. jlll . Q| . Spins were initially aligned along the z direction, 
in parallel with both the easy axes and applied fields. Stochastic simulation was performed 
on this initial configuration for 100 {'joHk)~^, in order to achieve the quasi-equilibrium state. 




FIG. 2: (color online) Dispersion relation for the simulated spin wave mode. Simulation parameters 
are: chain length N = 200, free boundary condition, thermal condition KyV/hBT = 50, exchange 
coupling strength J/2KuV = 1 and damping constant a = 0.1. Kittel's model refers to the 
theoretical dispersion relation of Eq. Q . 

Space and time Fourier transforms were then performed on the off-axis components. From 
the resulting spin wave spectra, the peak frequency uo determined for a range of wavevector 
k. The resulting dispersion relation in Fig. (j21) shows a very good convergence between 
the simulated results (calculated from both LLG and MC) and the theoretical prediction of 
Eq. (0). 

However, there do exist limitations to the TQMC method. The TQMC method cannot be 
easily extended for some advanced micromagnetic simulations, which include e.g. the spin 
torque effect, as compared to the stochastic LLG model. In addition, even though we found 
the TQMC algorithm to be typically 2-5 times faster than the LLG-based simulation with 
an equivalent time step, it is still too inefficient to model long-time magnetization reversal 
of up to, say, 1 second. Nevertheless, our analysis has raised the possibility of developing 
advanced time-quantiable Monte Carlo methods, based on e.g. the N-fold way Monte Carlo 
algorithm and kinetic Monte Carlo method j^, for micromagnetic studies. 



VI. CONCLUSION 



We have derived the time quantification factor which relate the time-scales of TQMC 
method and Langevin dynamics in the stochastic simulation of an interacting array of 
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nanoparticles. The time quantification factor is found to be the same as that derived pre- 
viously for isolated single-domain particles, up to the linear-order in the time-step size At. 
No explicit form of the Hamiltonian is implied in the derivation, which means that the 
equivalence between the two stochastic schemes is general and independent of many physi- 
cal and simulation parameters. To demonstrate this, we implement the TQMC scheme for 
the study of i) time evolution and magnetization reversal in a square spin array, and ii) 
spin wave behavior in a one-dimensional interacting array of particles. In the case of (i), 
we show a close correspondence between the TQMC and LLG results for a wide range of 
coupling strength. The equivalence remains valid even when the reversal mode changes as 
J is increased. In the case of (ii), the numerical verification of the time quantification factor 
is provided by the close agreement of the spin wave dispersion curves as obtained by both 
TQMC and Langevin dynamical methods. The two curves also show a very close agreement 
to the known theoretical dispersion relation. Our analytical derivation and numerical stud- 
ies thus justify the applicability of the TQMC method for stochastic micromagnetic studies 
in most cases where the LLG equation applies. 



VII. APPENDIX 



A. FP coefficients for the LLG equation 

A previous study of the effect of thermal fluctuations in an interacting spin array system 
based on the Langevin scheme, showed that the inter-particle interactions do not result in 
any correlations of the thermal fluctuations However, to the best of our knowledge, 

a detailed derivation of the FP coeffcients for an interacting particle system has not been 
presented. Hence we include, as an appendix, a derivation of FP coefficients for an interacting 
particle system. We extend Brown's derivation for the FP coefficients of isolated single 
domain particles to obtain the FP coefficients for the case of interacting particles. 

The thermal field h(t) representing the thermal fluctuations, according to Brown has 
the properties of a white noise, i.e. 

m)) = 0, (h^hp)) = 2D ■ 5„5,,5{t) (8) 

where i,j = {1,2,3} denote the Cartesian coordinates component {x, and p,q = 
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{1, . . . , A^} refer to the p^^ and q^^ spin in the hst. Hence, if: 

ft+At 



Kf ^ h^dt' (9) 

then 

{Kf) = 0, {KfK]) = 2D-S,,S,,At (10) 
Rewriting Eq. in spherical coordinates, we obtain 

Left side = -j- = t^-tt + ^ jr = ■ 9i + ■ sm9iLpi (11) 

at aOi at oipi at 

, 1 70 f dE f dE 

Right side = — TTTTT^ 9V ^ ~^ h a ■ Sj X x 



h'—- '^——\e ( '—-h'-^—\e (12) 

dOi ^ sin 9i dipi J y ddi sin 6i dipi J ^'^ 

in which the partial differential relationships such as^ = ^|^ + ^^ = §f + f^^j^e^ 
have been used. With the inclusion of the thermal fluctuation, additional terms will be added 
into the right side as (sj x h(t) + a ■ Sj x (sj x h.{t))). By considering the relation 

between Cartesian and spherical base vectors: 

* = sin 6 cos (/? ■ + cos 6 cos ■ eg — sin ■ 
j = sin 6* sin ■ + cos 6^ sin ■ + cos ip ■ 

k = cos 6 ■ Br — sin 6 ■ eg (13) 
and equating Eqs. (fTTj) and (fT^ . we thus obtain 2N simultaneous equations as: 



at ' smfc'j 

^ = + h'^:^H' (14) 

dt ^ mi6i sin^^i ^ ^ 

where 

+ + (15) 

if^i and are the contributions of h(t) to the generalized forces corresponding to 9i and 

{2KuV)^^Hg- = h\{t) cos 9 i cos (fi + h2{t) cos 9 i sin (fi — hl(t) sin 9 i 
{2KuV)-^H^^ = -h\{t)sin9isinipi + h^^it) sin 9i cos (pi (16) 
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Eq. ()14|) can be expressed directly in a general form as: 

i^ = m^)+i2^i{^^hi{t) (17) 
fc=i 

where x represents the set of 2N variable {xf } (here i = {1,2} denotes angular coordinates 
{6, (f} and p = {1, 2, ... , N} refers to the p^^ spin in the list). To evaluate the FP coefficients 
Ax^ and Bx^xj , we need Axj only to terms of the order At for A^. and only to terms of order 
(At)i/2 Q^^^^^ Taking note of Eq. ©, Axi itself is of order (At)^/^ Expanding F[{x) 
and G^i^{x) in Taylor's series at initial state Xq: 



Ffix) = Ff (xo) + Y.nf- M + I E ■ A^^A^F + • • • 
GU^) = Glixo) + E ■ ^ E ■ ^^I^^F + ■ • ■ (18) 

<},j q,r,j,i 

where, for example, F[j^ = dFf jdx^^ and G\'i!^ = dG^j^/dx''-. Hence by integration of Eq. (fTTj) 
with respect to At, and truncate the terms that has order higher than At, we have: 

Axf = Ff At + ^ J hl{t,)dt^ + J2GTj ^x]hl{U)dt, (19) 

and in the last integral we may express AxJ to the order of At^/^, namely, as 
i:,G],C'hj{t,)dh. Thus, 

A< = Ff At + Y.GI hl{h)dh + E ^^!;-^'^ / / Kiti)hKt2)dh (20) 

the second term is of order At^/^, the others of order At; therefore, to the first order in At: 

/■Ai /.At 

ax^ax] = j2gig'^J dh hi{t,m,)dh (21) 

k,i -^0 Jo 

It is easily seen that the double integral in Eq. (j20j) is half that in Eq. (j21|) . We now 
evaluate the statistical average by considering Eq. (|10|) and dividing by At: 



= lim = F[ + D-y G^ip% 

At-^O At ^ ^ 



k 



(Ax^Ax'^) 

= jim = ■ E ^'fc^'fc ■ ^P'^^^^ (22) 



10 



In the present application, 



,v udE ^, 1 BE 



dOp sin 9p d(pp 

pp = g'^^ - h'^—^ (23) 

^ sin 9p dOp sin^ 6p d(pp 



and 



{2KuV) ^Gii = h' cos 9p cos (fp — g' sin (fp 
{2KuV)^^Gi2 = h' cosOpSimpp + g' cosipp 

(2ir„\/)-iG?3 = -h'sinOp (24) 
{2KuV)~^G2i = —g' cot 9p cos ipp — h' esc Op sin Lpp 
{2KuV)~^G22 = cot 6'p sin + /i' CSC 6*^ cos 
(2ir„F)-iGf3 = g' 

Partial differentiation of Eqs. ()24|) with respect to Op and (pp gives the formulas for the twelve 
quantities G^^i^- = l,2;k = 1,2,3). Substitution of the values of F[, G^f^ and G^j^'j into 
Eqs. fl22|) gives the value of FP coefficients for LLG dynamical equation as follows: 

A^^^ = -h'- g' hA;'cot6'i 

oOi sin Oi otpi 

,LLG ^ , 1 dE 1 dE 

^ smO.dO, sm^Oidipi 

B^i^ = 2k' . 5,, (25) 



_ _J_2k' . Si, 
'^^'^^ sin^ Oi ^ 

tdLLG -ryLLG r\ 

where k' = D{h'^ + g''^){2KuVY is to be determined since the value of D is still unknown. 
Substituting Eqs. into Eq. (0) and taking note that P({6'}, {v^},^) should reduce to the 
Boltzmann distribution at statistical equilibrium {dP/dt = 0), one thus obtain the value of 
k': k' = h'/(3. 

B. FP coefficients for TQMC 

We next derive the FP Coefficients for TQMC. The Monte Carlo algorithm starts with a 
random selection of the spin site. We consider the z*^ spin in the list. For a trial move with 
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the displacement vector to be of size (r^ < R) and angle with respect to e^, we have 
the corresponding change with respect to 6i and ipi as 

A9i = —Ti cos ctj + cot 6i siv? ai + O (rf ) 

Aifi = Vi^^ sin ai + r,^ cos q;^ sin + O (rf ) (26) 



The displacement probability of the size to be is given by Nowak et al. |l| as 



pij-i) = 3^R^ - rf/27iR^ (27) 

and the acceptance probability for this trial move is given by the heat bath rate as 

1 



A{AE) 



1 + exp {f3AE) 
U^_l (BE.. BE 



2V 2^[W^'-^W.^^'^^ '''' 



where AE is the energy change in the random walk step and (3 = {kBT)^^. Integrating over 
the projected surfaces [see Fig. (1) in Ref. ^1 for a clear diagram], we obtain a series of the 
required mean 

/•2n i-R T>2 arp 

{A9,) = / daj {rdn)Ae,-p{n)-A{AE) = — {cote, -13— ) + 0{R') 
Jo Jo 2U aUi 

sm t/j 20 0(pi 

{A9^) = ^ + 0{R') (29) 

sm 6i 20 
{Ae,A^,) = 0{R') 

Let subscript i (j) refers to the i*^ (j**^) spin in the list and X, Y denote either 9 or ip. 
One easily finds that when i ^ j: (AXiAYj) = 0. This is because in the Monte Carlo 
algorithm, only 1 spin site is chosen at each Monte Carlo step. Truncating the higher order 
terms in the above equations and including the probability factor of (1/A^) in choosing the i^^ 
spin from all spins, we then obtain the FP coefficients for TQMC method as in Eqs. 
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